Aminet 33/Aminet 33 - October 1999.iso/Aminet/misc/math/TCalcStats2c.lha/TCalcStats2c/AREXX/Normality.rexx

OS/2 REXX Batch file | 1999-07-03 | 12.9 KB | 697 lines

/* Chi-Sq Goodness of Fit Normality Test */

options results
if ~show('P','TCALC') then do
    address command 'run turbocalc:turbocalc'
    address command 'waitforport TCALC'
    loadflag=1
end
address 'TCALC'
'DEFPUBSCREEN()'
/* Add-in Rexx Math Library needed for some routines */
signal on syntax
if ~show('l','rexxmathlib.library') then
   call addlib('rexxmathlib.library',0,-30)
if ~show('l','rexxreqtools.library') then
   call addlib('rexxreqtools.library',0,-30)
if ~show('l','rexxsupport.library') then
   call addlib('rexxsupport.library',0,-30)
  /* add to library list */
signal off syntax

/* Start Main Routine */
if loadflag=1 then 'Load()'
'ActivateWindow()'
range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
colon=pos(":",range)
if colon=0 then do
    'Message "Please select a range before executing this script"'
    'DEFPUBSCREEN("Workbench")'
    exit
end

/* Find cell references and cell, column numbers */
start_cell=substr(range,1,colon-1)
end_cell=substr(range,colon+1)
start_row=cellrow(start_cell)
end_row=cellrow(end_cell)
start_col=cellcol(start_cell)
end_col=cellcol(end_cell)
NRows=end_row-start_row+1
NCols=end_col-start_col+1
If NCols>1 then do
    'Message "Only one column of data allowed!"'
    'DEFPUBSCREEN("Workbench")'
    exit
end

/* Get cell reference for output range */
out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
if out_cell="" then do
    'DEFPUBSCREEN("Workbench")'
    exit
end
if length(out_cell)<2 | datatype(left(out_cell,1),'n')=1 then do
    'Message "Invalid cell reference"'
    'DEFPUBSCREEN("Workbench")'
    exit
end
BW=rtgetstring("0","Enter Desired Class Interval: 0 = Calculate","Input Request")
BW=strip(BW)
if rtresult=0 then DO
    'DEFPUBSCREEN("Workbench")'
    exit
end
min=rtgetstring("0","Enter First Lower Class Boundary: 0 = Calculate","Input Request")
min=strip(min)
if rtresult=0 then DO
    'DEFPUBSCREEN("Workbench")'
    exit
end
/* Suppress Screen Redraw to Speed Things Up */
'Refresh 0'

/* Open a small output window on tcalc screen*/
fo=0
CR='0a'x
DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
     call writeln(6Info, DisplayMsg)
    fo=1
end
else do
    'Message "TCALC Screen not available for Progress messages"'
end
CALL DELAY(150)

/* Get cell references for top cell in each column */
'SelectCell' start_cell
do col=start_col to end_col
    'GetCursorPos'
    top_cell.col=result
    'Column 1'
end

/* Get labels for later use on output */
'SelectCell' start_cell
'GetValue'
testlabel=result
testlabel=strip(testlabel)
if datatype(testlabel,'n')=1 then do
    labelflag=0
    title="Column "||x
end
else do
    labelflag=1
    NRows=NRows-1
    title=testlabel
end
if fo then call writech(6Info,"Progress...10 ")

/* Get data from cell range */
col=start_col
lav=0
tot=0
count=0
total=0
do x=1 to NCols
    'SelectCell' top_cell.col
    if labelflag=1 then 'CursorDown 1'
    do y=1 to NRows
        'GetValue'
        valtest=result
        if datatype(valtest)='NUM' then do
            'GetValue'
            val=result
            val=strip(val)
            data.y=val
            total=total+val
            count=1+count
        end
        'CursorDown 1'
    end
if fo then call writech(6Info,"20 ")

/* Calculate Mean and Standard Deviation */
mean=0
mean=total/count
dat=0
sum=0
do y =1 to count
    dat=data.y
    sum=sum+(dat-mean)**2
end
N=count-1
var=(sum)/N
sd=sqrt(var)
if fo then call writech(6Info,"30 ")

/* Sort Values */
call Sort()
if fo then call writech(6Info,"40 ")

/* Calculate Minimum, Maximum*/
max=0
N=count
if min=0 then min=trunc(data.1)
max=trunc(data.N+.5)
if fo then call writech(6Info,"50 ")

/* Calculate lcb and ucb Range */
cr='0a'x
test=1
x=1
Nbins=0
bin.=0
ucb.=0
if BW=0 then DO
    Nbins=trunc(sqrt(count))
    if Nbins<5 then Nbins=5
    if Nbins>15 then Nbins=15
    BW=(max-min)/Nbins
End
Else DO
    Nbins=trunc((max-min)/BW)
    if Nbins<5 then Nbins=5
    if Nbins>15 then Nbins=15
End
if BW<.5 then cf=.0001
if BW<1 & BW>=.5 then cf=.001
if BW<10 & BW>=1 then cf=.01
if BW<100 & BW>=10 then cf=.1
if BW>=100 then cf=1
lcb.1=min
do x=2 to Nbins
    z=x-1
    lcb.x=(lcb.z)+BW
    ucb.z=(lcb.x)-cf
end
if fo then call writech(6Info,"60 ")

/* Calculate Mean Deviation etc */
meandev.=0
stscore.=0
Do y=Nbins-1 to 1 by -1
    meandev.y=ucb.y-mean
    stscore.y=(meandev.y)/sd
End
if fo then call writech(6Info,"70 ")

PBZ.=0
PWZ.=0
PBZ.Nbins=1
Do y=Nbins-1 to 1 by -1
    x=y+1
    call NPROB(stscore.y)
    PBZ.y=P
    PWZ.x=(PBZ.x)-PBZ.y
End
if fo then call writech(6Info,"80 ")

PWZ.1=PBZ.1
Do y=1 to Nbins
    ExpFreq.y=PWZ.y*count
End

/* Calculate Observed Frequencies */
x=1
freq.=0
do y=1 to NRows
    z=x+1
    if x<Nbins then
    if (data.y)>=(lcb.x) & (data.y)<(lcb.z) then freq.x=(freq.x)+1
    else do
        y=y-1
        x=x+1
    end
    else freq.x=(freq.x)+1
end

/* Calculate Residuals and Chi-square */
R.=0
Chi=0
D.=0
DD=0
Do y=1 to Nbins
    DD=(freq.y)-(ExpFreq.y)
    If DD=0 then DO
        D.y=0
        R.y=0
    End
    Else DO
        D.y=(DD)**2
        R.y=DD/sqrt(ExpFreq.y)
    End
    Chi=Chi+((D.y)/ExpFreq.y)
end
if fo then call writech(6Info,"60 ")

/* Degrees of Freedom? */
w=rtgetlong(0,"Change df based on parameters estimated? - default=0","Value Required")
if rtresult=0 then DO
        'DEFPUBSCREEN("Workbench")'
        exit
End

df=Nbins-1-w
if df<1 then df=1

/* Calculate Probability */
PC=CHIPROB(Chi,df)
if fo then call writech(6Info,"80 ")

/* Calculate Critical */
Chicrit1=CHICRIT(.95,df)
if fo then call writech(6Info,"90 ")
Chicrit2=CHICRIT(.99,DF)
if fo then do
    call writeln(6Info,"100")
    call writeln(6Info,"Writing output to window...")
end

/*Output*/
'SelectCell' out_cell
'ColumnWidth 18'
'Put "Chi-Sq Goodness of Fit of Normal Frequency Distribution"'
'CursorDown 2'
title=""""||title||""""
'Put' title
'CursorDown 2'
'GetCursorPos'
st_cell=result
'Alignment 2'
'Put "l.c.b"'
'Column 1'
'Alignment 2'
'Put "u.c.b"'
'Column 1'
'CursorUp 1'
'Alignment 2'
'Put "Mean"'
'CursorDown 1'
'Alignment 2'
'Put "Deviation"'
'Column 1'
'CursorUp 1'
'Alignment 2'
'Put "Standard"'
'CursorDown 1'
'Alignment 2'
'Put "Score"'
'Column 1'
'CursorUp 1'
'Alignment 2'
'Put "Proportion"'
'CursorDown 1'
'Alignment 2'
'Put "Within"'
'Column 1'
'CursorUp 1'
'Alignment 2'
'Put "Expected"'
'CursorDown 1'
'Alignment 2'
'Put "Frequency"'
'Column 1'
'CursorUp 1'
'Alignment 2'
'Put "Observed"'
'CursorDown 1'
'Alignment 2'
'Put "Frequency"'
'SelectCell' st_cell
'CursorDown 1'
Do y=Nbins to 1 by -1
    'Put' format(lcb.y,,2)
    'CursorDown 1'
End
'SelectCell' st_cell
'Column 1'
'CursorDown 2'
Do y = Nbins-1 to 1 by -1
    'Put' format(ucb.y,,2)
    'CursorDown 1'
End
'SelectCell' st_cell
'Column 2'
'CursorDown 2'
Do y=Nbins-1 to 1 by -1
    'Put' format(meandev.y,,2)
    'CursorDown 1'
End
'SelectCell' st_cell
'Column 3'
'CursorDown 2'
Do y=Nbins-1 to 1 by -1
    'Put' format(stscore.y,,4)
    'CursorDown 1'
End
'SelectCell' st_cell
'Column 4'
'CursorDown 1'
Do y= Nbins to 1 by -1
    'Put' format(PWZ.y,,4)
    'CursorDown 1'
End
'SelectCell' st_cell
'Column 5'
'CursorDown 1'
Do y=Nbins to 1 by -1
    'Put' format(ExpFreq.y,,2)
    'CursorDown 1'
End
'SelectCell' st_cell
'Column 6'
'CursorDown 1'
Do y=Nbins to 1 by -1
    'Put' format(freq.y,,2)
    'CursorDown 1'
End
'SelectCell' st_cell
'CursorDown' Nbins+3
'Put' "Chi-Square:"
'CursorDown 1'
'Put "d.f.:"'
'CursorDown 1'
'Put "P(CHI<=chi):"'
'CursorDown 1'
'Put "Chi-Critical (95%):"'
'CursorDown 1'
'Put "Chi-Critical (99%):"'
'CursorUp 4'
'Column 1'
'Put' format(Chi,,4)
'CursorDown 1'
'Put' df
'CursorDown 1'
'Put' format(PC,,6)
'CursorDown 1'
'Put' format(ChiCrit1,,4)
'CursorDown 1'
'Put' format(ChiCrit2,,4)
'Refresh 1'
'Refresh 2'
/*'Message' "Finished"*/
/*indicate the main script is finished*/
DisplayMsg="Cleaning up ...."||CR||"Exiting"
result=writeln(6Info, DisplayMsg)
if result~=0 then do
    /*Wait 3 seconds*/
    CALL DELAY(150)
    /* close window*/
    result=close(6Info)
end
'DEFPUBSCREEN("Workbench")'
exit

/* Procedures */

cellrow: procedure
do
    parse arg cell
    do charpos=2 to length(cell)
    if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
    end
    return 0
end
Return

cellcol: procedure
do
    parse arg cell
    labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
    cell=upper(cell)
    len=length(cell)
    val=0
do charpos=1 to len
    if datatype(substr(cell,charpos,1),n) then
    do cell=reverse(substr(cell,1,charpos-1))
    do x=1 to length(cell)
    val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
    end
    return val
    end
    end
    return 0
end
Return
/* It is important to put the exposed array at the end of the next line */
Sort: procedure expose  NRows data.
L=(xtoy(2,int(log(NRows)/log(2))))-1
    Do Until L<1
    L=trunc(int(L/2))
    Do J=1 to L
            Do K=J+L To NRows By L
            I=K
            dumdat=data.I
            Do while I>L
                y=I-L
                If data.y ~> dumdat then Leave
                data.I=data.y
                I=I-L
            End
            data.I=dumdat
            End
        End
    End
Return

syntax:
     if arg(1)='FAIL' then do
        'Message "Library is unavailable."'
        'DEFPUBSCREEN("Workbench")'
        exit
        end
    'DEFPUBSCREEN("Workbench")'
    exit



Format:  procedure

     arg number, before, after
     CallLine = SIGL
     if ~datatype(CallLine, 'N') then CallLine = '??'

     /* Make sure we have a number as first (required) argument    */
     if ~datatype(number, 'N') then do
        if number = '' then
           fc = 17     /* Wrong number of arguments           */
        else
           fc = 47     /* Arithmetic conversion error             */
        signal FormatSyntaxError
     end
     num = number + 0
     if before = '' & after = '' then
        return num
     else do
        parse var num integer '.' fraction
        if before = '' then before = length(integer)
        if after = '' then after = length(fraction)
        if ~datatype(before, N) | ~datatype(after, N) then
           do fc = 18
           signal FormatSyntaxError
       end
        if before < length(integer) then do
           fc = 18
           signal FormatSyntaxError
        end
        if after ~= length(fraction) then do
           fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
        if integer<1&integer>-1 then integer=integer
           else integer = integer + (fraction % 1)
           fraction = substr(fraction, 3)
        end
        if fraction >= 0 then
           return right(integer, before)'.'fraction
        else
           return right(integer, before)
     end

 FormatSyntaxError:
        if show('F', STDERR) then
           call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
        else
           mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
        'Message' mess
        parse source Func .
        if Func = 'FUNCTION' then do
        'DEFPUBSCREEN("Workbench")'
           exit "Err"
        end
        else do
        'DEFPUBSCREEN("Workbench")'
           exit 10
        end

NPROB: Procedure Expose P Q PDF

    arg Z
    A0=0.5
    A1=0.398942280444
    A2=0.399903438504
    A3=5.75885480458
    A4=29.8213557808
    A5=2.62433121679
    A6=48.6959930692
    A7=5.92885724438
    B0=0.398942280385
    B1=3.8052E-8
    B2=1.00000615302
    B3=3.98064794E-4
    B4=1.98615381364
    B5=0.151679116635
    B6=5.29330324926
    B7=4.8385912808
    B8=15.1508972451
    B9=0.742380924027
    B10=30.789933034
    B11=3.99019417011
    ZABS = ABS(Z)
    Y = A0*Z*Z
    PDF = EXP(-Y)*B0
SELECT
    WHEN (Z>=-1.28) & (Z<=1.28) then DO /*Z BETWEEN -1.28 AND +1.28*/
        Q = A0-ZABS*(A1-A2*Y/(Y+A3-A4/(Y+A5+A6/(Y+A7))))
        IF (Z<0) then do
            P = Q
            Q = 1.0-P
        End
        Else P = 1.0-Q
        RETURN
    End
    WHEN (Z>1.28) & (Z<=12.7) then do /*ZABS BETWEEN 1.28 AND 12.7 */
        Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/(ZABS+B9+B10/(ZABS+B11))))))
        IF(Z<0) then Do
            P = Q
            Q = 1.0-P
        End
        Else P = 1.0-Q
        RETURN
    End
    WHEN (Z<-1.28) & (Z>=-12.7) then do /*ZABS BETWEEN 1.28 AND 12.7 */
        Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/(ZABS+B9+B10/(ZABS+B11))))))
        IF(Z<0) then Do
            P = Q
            Q = 1.0-P
        End
        Else P = 1.0-Q
        RETURN
    End
    OTHERWISE /*Z FAR OUT IN TAIL*/
        Q = 0
        PDF = 0
        IF(Z<0) then do
            P = Q
            Q = 1.0-P
        End
        Else P = 1.0
        RETURN
End
Return

CHIPROB: PROCEDURE

    PARSE ARG X2,K1
    X3=0.5*X2
    K3=0.5*K1
    X=K3+1
    LG=LOGGAMMA(X)
    S=0
    A=EXP(K3*LN(X3)-LG-X3)
    IF A~=0 THEN DO
        T=1
        S=1
        G=K3
        DO WHILE (T/S) >.000001
            G=G+1
            T=T*X3/G
            S=S+T
        END
    END
    PO=S*A
RETURN PO

CHICRIT: PROCEDURE

    PARSE ARG PO,K1
    AO=.0001
    E=.005
    E2=E+E
    PA=PO
    SELECT
    WHEN K1=1 THEN DO
        PO=.5+PA/2
        Z=NORMALPP(PO)
        X1=Z*Z
    END
    WHEN K1=2 THEN DO
        X1=1-PA
        X1=-2*LN(X1)
    END
    OTHERWISE
        Z=NORMALPP(PO)
        A1=2/(9*K1)
        W=1-A1+Z*SQRT(A1)
        XO=K1*(W**3)
        DO FOREVER
            X2=XO
            PO=CHIPROB(X2,K1)
            F1=PO-PA
            X2=XO+E
            PO=CHIPROB(X2,K1)
            F2=PO
            X2=XO-E
            PO=CHIPROB(X2,K1)
            F2=F2-PO
            F2=F2/E2
            X1=XO-F1/F2
            IF ABS(X1-XO)>AO THEN XO=X1
            ELSE LEAVE
        END
    END
    X2=X1
RETURN X2

LOGGAMMA: PROCEDURE

    ARG XA
    C1=76.18009173
    C2=-86.50532033
    C3=24.01409822
    C4=-1.231739516
    C5=.001208580
    C6=-.000005364
    C7=2.506628275
    X1=XA-1
    W=X1+5.5
    W=(X1+.5)*LN(W)-W
    S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
    S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
    L=W+LN(C7*S)
RETURN L

NORMALPP: PROCEDURE

    ARG P0
    A1=2.515517
    A2=.802853
    A3=.010328
    A4=1.432788
    A5=.189269
    A6=.001308
    Q=.5-ABS(P0-.5)
    SN=SIGN(-2*LN(Q))
    IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
    ELSE W=SQRT(-2*LN(Q))
    W1=A1+W*(A2+A3*W)
    W2=1+W*(A4+W*(A5+A6*W))
    ZZ=W-W1/W2
    SELECT
        WHEN (P0-.5)<0 THEN TT=-1
        WHEN (P0-.5)=0 THEN TT=0
        otherwise TT=1
    END
    ZZ=ZZ*TT
RETURN ZZ